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Miscellaneous magnetic systems are being recently intensively investigated because of their po- 
tential applications in modern technologies. Nonetheless, a many body dynamical description of 
complex magnetic systems may be cumbersome, especially when the system exhibits a geometrical 
frustration. This paper deals with simulations of the classical XY model on a three dimensional 
triangular lattice with anisotropic couplings, including an analysis of the phase diagram and a Bo- 
goliubov description of the dynamical stability of mean-field stationary solutions. We also discuss 
the possibilities of the realization of Bose-Hubbard models with complex tunneling amplitudes in 
shaken optical lattices without breaking the generalized time-reversal symmetry and the opposite, 
i.e. real tunneling amplitudes in systems with the time-reversal symmetry broken. 
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I. INTRODUCTION 

This paper was inspired by an experiment of the Ham- 
burg University group who managed to build the first ex- 
perimental realization of the two-dimensional (2D) quan- 
tum simulator of classical frustrated magnetism in an 
ultra-cold atomic gas in a triangular optical lattice 0- 
3]- Frustrated magnetism is the result of competition 
between interactions and geometry of a lattice and con- 
stitutes a very active field of research 0-0] • 

Quantum simulators are easily controllable quantum 
systems that can be employed to mimic others Q. If 
considered systems are similar enough, one may be able 
to map one system's Hamiltonian onto the Hamiltonian 
of the simulator. Consequently, a time evolution of the 
quantum simulator imitates a simulated one. Quantum 
simulators are extremely useful when issues in consider- 
ation are substantially too troublesome for computers or 
their direct observation is laborious or hardly possible. 

An exquisite example of a quantum simulator is a sys- 
tem of ultra-cold bosonic gas in an optical lattice |8rlll| . 
Among others, that is because of a great flexibility of 
optical lattices, as their geometry can be changed easily 
and a possibility to manipulate atomic interactions by 
changing depth of a lattice or via Feshbach resonances 
[l2l fig ). An optical lattice has a structure of an ideal 
crystal, being devoid of defects present in material crys- 
tals, and therefore is a natural simulator of many solid 
state physics problems, such as miscellaneous spin mod- 
els. 

In our work, we simulate the classical XY model with 
nearest neighbor interactions characterized by Jy cou- 
plings 



H 



where 



S Ti = (cos 6 r i , sin 9 Ti ) , 



(1) 



(2) 



is a two component classical vector representing a direc- 
tion of a spin on i-th lattice site. Simulations of such 
a model in ultra-cold atomic gases are possible due to 
an identification of the directional angle 9 Ti with a local 
phase of bosonic condensate wavefunction 9 ri = k ■ of 
atoms in an optical lattice of the same geometry as the 
spin lattice. The vector k is a quasi-momentum vector 
of a single atom and the couplings Jij are related to el- 
ements of tunneling amplitude matrix. By doing so, an 
energy per a single spin can be associated with a dis- 
persion relation of an atom of an ideal non-interacting 
gas. A key to simulations is a possibility of manipulating 
Jij values. One may conclude that it is impossible to 
simulate anti-fcrromagnctism since tunneling amplitudes 
Jij naturally have to be non-negative. However, it turns 
out that they can be effectively made negative or even 
complex 



they can be 
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A model in our interest is the classical XY model on a 
3D triangular lattice corresponding to the trigonal struc- 
ture [Tijj with basis vectors forming equal angles of it/ 3 
and anisotropic tunneling amplitudes. The primitive vec- 
tors of the considered Bravais lattice read 



ai = ae x , 



a 2 
a.3 
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where a is the lattice constant. The vectors ai,a2,a3 
form a regular tetrahedron, see Fig. [1] 

In the present publication we mainly focus on real tun- 
neling amplitudes which have three different values de- 
pending on a tunneling direction as depicted in Fig. [T] 
Following the preceding identification, the energy per 
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II. CONSTRUCTION OF 3D TRIANGULAR 
LATTICE MODEL 




Cli 



Figure 1: (Color on line) Triangular optical lattice corre- 
sponding to the trigonal Bravais lattice with anisotropic tun- 
neling amplitudes. Lattice vectors a; form a regular tetra- 
hedron. Amplitudes Ji related to tunneling along different 
directions are indicated in the figure, i.e. J\ corresponds to 
the tunneling along dashed (red) lines, J2 along solid (blue) 
lines, and J3 along long-dashed (green) ones. 



spin depends only on one vector k and adopts a form 



E(\s) = -2| Ji cos (k • ai) + J 2 cos(k-a 2 ) + 

+ cos(k ■ (a 2 — ai)) + J 3 cos (k ■ a 3 ) + (4) 
+ cos (k ■ (a 3 - ai)) + cos (k • (a 3 - a 2 )) 



This paper is organized as follows. In Scc.[TJJwe present 
a background on how to make simulations of classical 
frustrated magnetism possible. Firstly, we propose a way 
of constructing a 3D triangular lattice. Then, to manip- 
ulate Jij we revise a procedure proposed in Refs. 0, [lj] 
to apply it to our model. In Sec. IIIII we describe the re- 
sults of our simulations, i.e. a phase diagram and phase 
transitions. In Sec. IIVI we analyze the stability of sta- 
tionary states within the Bogoliubov formalism. Sec- 
tion |V] is devoted to a discussion of the relation between 
complex tunneling amplitudes in effective Bosc-Hubbard 
models and the presence or absence of the generalized 
time-reversal symmetry of systems. Especially we show 
that complex amplitudes can be realized by the shaking 
of an optical lattice that does not break the generalized 
time-reversal symmetry. Finally, in Sec. IVII we conclude. 



A. Experimental setup 

Optical lattices used to be formed by a system of uni- 
polarized counter-propagating laser beams that form a 
standing-wave. As practical as it is for quasi-cubical lat- 
tices, it is hardly useful when dealing with triangular 
ones. The interfering running waves can be shifted to 
coincide at any angle. In the case of 2D triangular lat- 
tice three coplanar rcd-dctuncd and uni-polarizcd beams 
should be pointed inwards forming an angle of 2-7r/3 with 
each other. Such a construction was successfully real- 
ized experimentally in Refs. [H, This method could 
be further generalized into three-dimensions using four 
beams, however there emerges a problem of finding their 
proper polarizations. Under this circumstances, in or- 
der to construct a 3D triangular lattice, we propose to 
combine both methods. Namely, we first construct a 
2D lattice, which is in fact a lattice of one-dimensional 
tubes 0, l20j . Then, we add an extra oblique pair of 
counter-propagating beams that form a standing-wave. 
We choose, in a Cartesian basis, wave vectors of running 
waves (vectors ki,k 2 ,K3) and beams forming a standing 
wave (vectors ±k) to be 




n 2 = 




(5) 



where kl, k l are moduli of the wave vectors and n' L = 
kl/V3- The coordinate frame was chosen accordingly to 
construct the lattice given by the primitive vectors ([3]), 
where we imposed a condition a = 3/(V2kl) relating 
Kl and the lattice constant a. These proceedings are 
illustrated in Fig. fj] 



B. Effective description of driven lattices 

As it has been already sketched, the simulations are 
based on manipulations of tunneling amplitudes. It can 
be done without changing the lattice constant. An in- 
genious idea, proposed by Eckardt, Weiss, and Holthaus 
[14| , is to manipulate the values of tunneling matrix ele- 
ments by introducing a periodic lattice modulation. Ex- 
perimentally the motion of the lattice can be induced by 
a periodic frequency difference between any two interfer- 
ing beams. In result, one can make the lattice follow any 
desired trajectory R(f). 

An ultra-cold bosonic gas in a driven lattice is de- 
scribed in the co-moving frame by an explicitly time de- 
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Figure 2: (Color on line) Left panel illustrates the idea of real- 
ization of a 3D triangular optical lattice by cutting a 2D trian- 
gular lattice consisting of ID tubes. Arrows indicate primitive 
vectors of the resulting 3D lattice. Right panel shows config- 
uration of laser beams that create a 3D triangular optical 
lattice. Wave vectors «, and k of the laser beams correspond 
to Eq. ©. 



pendent Hamiltonian of the Bosc-Hubbard form 

-J2n l r l -F(t), (6) 

i 

where U is on site interaction energy and F(t) = -mR(t) 
is an inertial force [l4[ • We assume here that the shaking 
of the lattice does not lead us out of the lowest energy 
band of a lattice problem. It is valid provided that energy 
scale related to the driving frequency u> is considerably 
lower than the energy gap between the bands [2l| . 

Although energy eigenvalues do not exist for time de- 
pendent systems, an ultra-cold bosonic gas in a periodi- 
cally driven lattice does not constitute a non-equilibrium 
system. It was proven that every solution of the pe- 
riodically time dependent Schrodinger equation can be 
unambiguously expressed by the so-called Floquet states 
[2ll - |24j . That is the reason why the periodicity is of vital 
importance here. 

The term ^i^ 1 "* ' can be eliminated from the 
Hamiltonian ([5]) with the help of the unitary transforma- 
tion 

U(t)=eiq>(-^hiW i (t)J, (7) 




Figure 3: (Color on line) Panel (a) shows a configuration of 
spins in a 3D triangular lattice corresponding to the ground 
state of the system with tunneling amplitudes: Ji = — 0.3J, 
J2 = 0.2J and J3 = 0.3J. It belongs to a spiral phase indi- 
cated in Fig. [4] as S2. By rotating spin vectors ([2]) in every 
second aia2 plane by n we obtain the configuration shown in 
panel (b) which corresponds to the ground state of the system 
where Jj, — > — J3. The change of the sign of the J3 coupling 
changes the corresponding interactions from ferromagnetic to 
anti-ferromagnetic what, in turn, requires the rotation of the 
spins in order to obtain a ground state configuration. 



amplitudes 

J' l3 (t) = J l0 cxp U W^t)) , Wij(t) = Wi(t) - Wj(t). 

(9) 

If we assume now a high frequency regime, i.e. when 
fku is considerably larger than all energy scales in the 
Hamiltonian, then we can time average fast oscillating 
phases and our system can be described by the effective 
time independent Bose- Hubbard Hamiltonian 0, [2l[ 

6 = -e-^^+Ie^- 1 )- ( io ) 

(id) • 

with effective, renormalized tunneling amplitudes 

J!f = ^ £ dt J^t), (11) 

where T = 2ir/ui. Such effective models have been real- 
ized in several experiments 0, H3, • 



with 



Wi(t) = -n- f dt' F(f ). (8) 
Jo 



C. Effective tunneling amplitudes for 3D triangular 
lattice 



We consider a 3D triangular lattice corresponding to 
The resulting Hamiltonian has time dependent tunneling the trigonal Bravais lattice with primitive vectors forming 
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a regular tetrahedron that can be created by means of 
laser beams. 

In the absence of particle interactions (U = 0), en- 
ergy eigenvalues of the Hamiltonian (|10[) correspond to 
the dispersion relation of a single atom in the presence 
of a lattice potential. In the presence of the repulsive 
interactions (U > 0) and for a homogeneous system (i.e. 
(fii) = n where n is the mean number of atoms per lattice 
site), the ground state of the Gross-Pitaevskii equation 
[28| is related to the minimum of the dispersion relation. 
Now we would like to show how to realize cold atom ex- 
periments where the dispersion relation of a single par- 
ticle in the 3D triangular lattice is given by Eq. Q and 
the tunneling amplitudes J\ , J2 and J3 can be controlled 
exp er iment ally. 

When J3 coupling vanishes, the model is equivalent 
to two-dimensional one, which was already investigated 
[l|-[3|]. In the 2D case, the shaking of the lattice on an 
elliptical trajectory was used to manipulate the tunnel- 
ing amplitudes J\ and J2 0. ^ t urns out that in the 
3D case, changing the orientation of the ellipse and its 
shape allows us also to manipulate all three amplitudes. 
We assume that the triangular lattice is moving on an 
elliptical trajectory with the velocity 



R(t) 



[F x e x cos(u>t) - (F y e y 



F z e z )sm(ujt)}e{t), 
(12) 



where the Hcaviside step function 9(t) indicates that the 
shaking of the lattice begins at t — 0. Consequently, the 
inertial force F(£) = —mR(t) contains Dirac-dclta con- 



tribution which has to be taken into account [21| . In 
the case of the 3D triangular lattice under consideration, 
if we want to have only three different tunneling am- 
plitudes as indicated in Fig. [TJ the z-component of the 
inertial force has to fulfil 



p2 _ F 2 
p _ v x 

Then, the three tunneling amplitudes read 



(13) 



Ji 
■h 



J Jo 



F x a 



JJ °(^ F * +3F y 2 )> 

J Jo 



F 2 + F 2 



F u 



(14) 



where J is the tunneling rate of atoms between neigh- 
boring sites of the static (not shaken) triangular optical 
lattice and Jo represents the Bessel function of the zero 
order. The oscillatory character of the Bessel function 
allows us to change values and signs of the tunneling am- 
plitudes by changing the parameters of the shaking, i.e. 
F x , F y and a/ui. 

In Eq. (fl2|) we have chosen a particular direction of 
the velocity vector at the moment we turn the shaking 



on. By changing uit — >• u)t + ip, one can choose a different 
initial direction. However, we would like to stress that 
the phase ip is irrelevant and the tunneling amplitudes in 
Eq. (TTJJ are the same independently of the choice of ip. 

We have applied the unitary transformation (JT)) where 
the inertial force corresponds to the elliptical trajectory 
defined in Eq. (|12p and the resulting tunneling ampli- 
tudes (fT4"f are real valued. If we apply slightly different 
unitary transformation, i.e. with Wi{t) — > Wi(t) + 7.;, 
where 74 's are arbitrary real numbers, the resulting ef- 
fective Hamiltonian has also a form of the Bosc-Hubbard 
Hamiltonian (|10p but tunneling amplitudes may become 
complex valued. The freedom in choice of 7, (the so- 
called gauge freedom [T3|) is related to the freedom in 
choice of the global phase of the Wannier function local- 
ized at that site. Eigenstates of the system are obviously 
the same regardless of the choice of 74 's — for different 
sets of 7, 's they are written in different basis only. In ex- 
periments momentum distributions of atoms correspond- 
ing to ground states are measured. Momenta of atoms are 
detected in the laboratory frame. The Hamiltonian © is 
related to the frame co-moving with an optical lattice but 
the unitary transformation ([7]), where F(t) = -mR(t) 
and R(f ) is given in Eq. (fT2j) , leads to the frame where the 
momentum distributions coincide with those measured in 
the laboratory frame in experiments plj . Therefore, the 
choice of the same ji for all lattice sites (e.g. 7$ = 0) 
that results in the Bosc-Hubbard Hamiltonian with the 
effective tunneling amplitudes (|14[) is the unique choice 
justified by the measurement procedure. 



III. THE RESULTS OF SIMULATION 
A. Phase diagram and phase transitions 

In this subsection we present a description of a phase 
diagram and phase transitions of the model with real cou- 
pling constants, expressed by Eq. (| 14[) . In our study not 
only did we obtain ground state configurations which can 
be considered as generalizations of 2D ones but also new 
magnetic phases. In addition, we observe some interest- 
ing phenomena at phase transition boundaries. We will 
skip a detailed description of ground statephases of the 
system that are similar to the 2D case [lHS] and rather 
concentrate on new phenomena that are present in the 
3D lattice only. 

A ground state solution of the Gross-Pitaevskii equa- 
tion of a homogeneous atomic gas in a 3D triangular lat- 
tice corresponds to a minimum of the dispersion relation 
(|4]). The appearance of the J3 coupling implies some 
novelties as compared to the 2D case. First of all, ev- 
ery ground state configuration in the 2D case has two 
3D counterparts. The latter can be transformed into one 
another by rotating spins ([2]) in every second &1&2 plane 
by 7r — Fig. [3] presents an example. It is due to a sym- 
metry of the dispersion relation (U) which remains un- 
changed under a simultaneous replacement J3 — > — J3 



5 



h/J= 0. 




(a) ik y a/2K (b) »*,a/2* (c) tk y a/27T 

it . i? • J 



J 3 /J=0A 




Figure 4: (Color on line) Phase diagrams for two differ- 
ent positive values of J3, i.e. J3 = 0.1J (top panel) and 
Jj, = 0.4 J (bottom panel). Parameter regions of ferromag- 
netic and rhomboidal phases are denoted by F and R, respec- 
tively. Between these regions there is an intermediate phase 
denoted by I. There are three regions of different spiral phases: 
S2, Sl a and SI' 3 . Transitions between all phases are continues 
except the transitions between Sl a and Si' 3 and between S2 
and I. 



and k — > k + ba/2 where b3 is a primitive vector in the 
reciprocal space that is orthogonal to the sl\Sl2 planes. 
Because of this symmetry, in discussing the phase dia- 
gram it is sufficient to restrict ourselves to positive values 
of the J3 coupling only. Phase diagrams for two different 
positive values of J3 are presented in Fig. |U 

In the 2D case by changing J2 we can switch from fer- 
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Figure 5: (Color on line) Projection (along z-axis) of mo- 
mentum distributions of atoms prepared in ground states of 
a 3D triangular optical lattice. Panels (a) and (c) are related 
to the ferromagnetic and rhomboidal phases, respectively. 
Panel (b) shows two degenerate intermediate states (one is 
depicted by full symbols and the other by open ones) which 
allow us to pass continuously from the ferromagnetic phase to 
the rhomboidal one. Tunneling amplitudes are the following 
Ji = 0.2J, J 3 = 0.3J and J 2 = -0.1J (a), J 2 = -0.11J (b) 
and Ji — —0.3 J (c). 



romagnetic (F) to rhomboidal (R) phases and this tran- 
sition is the first order transition [l|, Q ■ In the 3D case 
the transition is no longer discontinuous because a new 
intermediate (I) phase emerges. Spins can now reorga- 
nize smoothly from the ferromagnetic configuration to 
the rhomboidal one as illustrated in Fig. [5] It is a pure 
consequence of the geometrical frustration of the system. 
Namely, a competition between J3 and J2 couplings en- 
ables intermediate configurations. 

Back in two-dimensions, spiral configurations SI and 
S2 are the only ones that exhibit signs of frustrations 
[H, Q- In 2D there is a continuous transition between 
these phases and they are separated by a boundary phase 
of anti-ferromagnetic chains. Naturally, this boundary 
appears when J2 vanishes while J\ stays negative, i.e. 
where anti-ferromagnetic chains can be formed. In the 
3D system, because of the extra interaction, the corre- 
sponding boundary is shifted. Anti-ferromagnetic chains 
do emerge in 3D as well, but the chains are not indepen- 
dent because there are correlations between spins belong- 
ing to different chains. 

What is curious in 3D, the region of the SI spiral 
configuration of the 2D model is now divided into two 
regimes that we denote by SI" and Sl^. The transition 
between these new spiral phases is of the first order (even 
for infinitesimal J3 values) and it goes through a bound- 
ary configuration with non-trivial degeneracy which is 
illustrated in Fig. [6] When we start with the SI" phase, 
change the parameters and reach the boundary we end 
with a different ground state than when we do the same 
but start with the Sl^ phase, see Fig. [5] A Similar phe- 
nomenon takes place on the boundary between the S2 
phase and the intermediate phase I, i.e. this transition is 
also discontinuous. 

All phases presented in Fig. 2] are doubly degenerated 
except F and R. The reason for the degeneracy is the in- 
variance of the dispersion relation (j4|) under reflection 
in the quasi-momcntum space. That is, by changing 
k — > — k in a given ground state one obtains another 
ground state. In the case of the F and R phases such 
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Figure 6: (Color on line) Curves form a set of quasi-momenta 
k corresponding to the ground state energy of the system at 
the border between Sl a and SI phases, i.e. for Ji = J2 ~ 
— 0.15J and J3 = 0.1J. The transition between the phases is 
discontinuous. Starting from the Sl a phase and passing the 
boundary, the ground state of the system jumps from one of 
the points indicated by the red (dash) vectors to one of the 
points indicated by the blue (solid) vectors, which represent 
two degenerate ground states of the Sl^ phase. Note that 
in the figure we do not use the Cartesian coordinate frame 
but the frame corresponding to the primitive vectors of the 
reciprocal space. 



a transformation does not produces a new ground state 
because it reduces to a translation by a primitive vector 
of the reciprocal space and thus leads to the same state, 
see Fig. [5] 



B. Bose-Hubbard model with complex tunneling 
amplitudes 

We have discussed the Bose-Hubbard model (jTUJ) 
with effective time-independent tunneling amplitudes 
Jff which are real valued. However, it is also possible to 
realize complex Jfj. Usually it is done by applying the 
shaking of an optical lattice that breaks the time-reversal 
symmetry [l5|, [TtJ (see also [l|| for a different method) 
but it turns out that one can realize complex JfJ 1 even 
in the presence of this symmetry — see Sec. |V] 

The introduction of complex coefficients modifies the 
dispersion relation ^ and shifts its extrema. This mod- 
ification reduces to a replacement of all couplings Jj 
in with their absolute values \Ji\ and shifting ar- 
guments of the corresponding cosines by their complex 
phase <fi, i.e. Ji cos (. . .) — > \ Ji\ cos (. . . — ifi). While ^3 



can be trivially eliminated with the help of a translation 
k — > k — f^b 3 (b 3 is a primitive vector in the recipro- 
cal space), the remaining phases tp\, ip2 change structure 
of the dispersion relation essentially. In particular, the 
quasi-momentum-reversal symmetry, whose presence is 
the reason for the double degeneracy of the spiral and 
intermediate phases discussed in Sec. IIII A| is broken. 
Furthermore, complex couplings make it possible to take 
a detour around any phase boundary and consequently 
avoid any first order phase transition present in the real 
coupling case. 

IV. STABILITY OF STATIONARY STATES: 
BOGOLIUBOV ANALYSIS 

A mean field analysis is a practical and convenient tool 
when describing weakly interacting cold bosonic gases 
in optical lattices. Stationary solutions of the Gross- 
Pitaevskii equation of a homogeneous system described 
by the Bose-Hubbard Hamiltonian are given in form of 
Bio eh waves 

Mn) = e lk °- r % (15) 

where ko is a quasi-momcntum vector and N s stands for 
a number of lattice sites. Even though we assume weak 
repulsive atomic interactions (i.e. nil <C J with mean 
number of n atoms per lattice site) they may contribute 
to the dynamical instability of stationary states. For this 
reason, we perform stability analysis of mean-field solu- 
tions. 

The Gross-Pitaevskii equation corresponding to the 
Bose-Hubbard model ((T0|) can be obtained by switch- 
ing from the operators to c- numbers, i.e. hi — > hi, in 
the Hcisenbcrg equation for hi [28| . Linearization of 
this equation leads to standard Bogoliubov-de Gennes 
equations 28]. In the present situation in order to find 
stationary Bogoliubov modes we employ the following 
ansatz 



"k(r l ) 


gik-r, 


' U k e ,;ko r ' " 






Vk e- ik °- r< 



(16) 



The resulting Bogoliubov-de Gennes eigenvalue problem 
is represented by a block diagonal matrix in which diag- 
onal elements are 2x2 matrices. Diagonalization of the 
blocks gives the eigenvalues 

A± = 1 (A£ k - A£ k ± <5A k ) , (17) 

where 

£A k = yf (A£ k + A£ k )(4n[/ + A£ k + A£k), 
A£ k - £(k +k)-£(k ), 

A^ k - £(k -k)-£(ko), (18) 
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Figure 7: (Color on line) Lyapunov exponent, i.e. imaginary 
part of Ak, see Eq. (|19[1 . in the units of bare tunneling ampli- 
tude J as a function of k x and k y for the ferromagnetic state 
and for tunneling amplitudes in the regime of the S1 Q phase. 
The third quasi-momentum coordinate k z corresponds to the 
maximum of Im(Ak). Tunneling amplitudes are Ji = —0.2 J, 
J2 = — 0.1J and J3 = 0.2 J, while top (bottom) panel is re- 
lated to nU = 50 J (nil = 0.2 J). 



and .E(k) is the dispersion relation When for a cer- 
tain k vector, a quantity ((5Ak) 2 is negative, then the 
corresponding eigenvalues are complex. In such a case, a 
stationary state (fT5")) is not dynamically stable and tends 
to collapse into a dominant Bogoliubov mode, i.e. a mode 
related to an eigenvalue with the largest imaginary part. 

The stability of excited states may be a relevant ex- 
perimental problem, when one prepares the system in a 
ground state i/)q and then changes the tunneling ampli- 
tudes. If ipo is not the ground state of the new Hamilto- 
nian, it may collapse on time scale depending how strong 
the repulsive particle interactions are. 

Usually, it is the most convenient to prepare the trian- 
gular lattice system in the ferromagnetic state, i.e. in the 
Bloch wave (fT5|) with ko = Q. Assume that after the 



preparation, experimentalists change abruptly the tun- 
neling amplitudes switching to a regime of the spiral S1 Q 
phase. If the initial state is the ferromagnetic one, the 
eigenvalue expression (|17[) simplifies to A^ = ±Ak with 

A k = v/A£ k (2nUTAEkj, (19) 

where AE'k = E(k) — E(0) is negative for some k. For 
U = the ferromagnetic state is still a stable stationary 
solution of the Gross-Pitaevskii equation. On the other 
hand, if nU is larger than max(— AEfc/2), a small per- 
turbation of the ferromagnetic wave-function will grow 
exponentially in time and occupations of the Bogoliubov 
modes with k corresponding to the Sl a ground states 
increase the fastest. Figure [7] shows imaginary part of 
Ak (the so-called Lyapunov exponent) where we can see 
that the maximal Lyapunov exponent is related to two 
quasi-momentum vectors. These vectors correspond to 
two degenerate ground states of the S1 Q phase. We may 
expect that in an experiment one of the ground states of 
the Sl a phase will be predominantly populated — which 
one depends on an initial perturbation and may be dif- 
ferent in different experimental realizations. Signatures 
of such a spontaneous symmetry breaking have been ob- 
served experimentally in the 2D case @]. 

In a very weak interaction regime, namely when nU 
is smaller than max(— A£k/2), a situation is different. 
Dominant modes can be found from the condition 

AE k + nU = 0. (20) 

An example of such a situation is presented in Fig. [7] 
where we can see a family of the most unstable modes. 
Interestingly the Bogoliubov modes with k corresponding 
to the Sl a ground states are stable. Thus, there is no 
decay of the ferromagnetic state towards the new ground 
states. 

The analyzed effects may by observed in experiments. 
For instance, we take realistic conditions nU = 50 J with 
J = h x 6.5 Hz and calculate the stability of the fer- 
romagnetic state in the region of the Sl a phase. For 
Ji = -0.2J, J 2 = -0.1J and J 3 = 0.2J, the en- 
ergy difference max(— AEk/2) = 0.45J and as a re- 
sult the ferromagnetic state should predominantly col- 
lapse towards one of the ground states on time scale 
of 1/ (2 1 Ak I ) ps 11ms. On the other hand, if we con- 
sider smaller mean number of atoms per lattice site, e.g. 
nU = 0.2 J, we may expect the ferromagnetic state to de- 
cay predominately towards a multiple modes correspond- 
ing to the maximal Lyapunov exponent on time scales of 
390ms. 



V. COMPLEX TUNNELING AMPLITUDES 
AND TIME-REVERSAL SYMMETRY 

In Sec. [IT] we have considered an optical lattice driven 
by a periodic external force and obtained an effec- 
tive time-independent Bosc-Hubbard Hamiltonian (fTU|) 
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Figure 8: (Color on line) Bravais lattice points (black circles) 
and amplitudes Ji corresponding to tunneling from a lattice 
point to the nearest neighbors. 
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Figure 9: An example of a velocity of a shaken ID optical 
lattice with no time reversal symmetry that leads to an effec- 
tive Bose-Hubbard model with real tunneling amplitudes, see 
Eq. JHJ(. 



with real tunneling amplitudes (|14p. The original time- 
dependent Hamiltonian is invariant under a generalized 
time-reversal symmetry. In other words, there exists an 
anti-unitary operator 7" which commutes with the Hamil- 
tonian. In such a case, it can be easily shown how to 
construct a basis in which Hamiltonian is represented by 
a real symmetric matrix (29j . As a consequence, a stan- 
dard idea to obtain effective Bose-Hubbard models with 
complex tunneling amplitudes is to apply external driv- 
ing that breaks that symmetry (l5l 1171. Hq. However, we 
will show that an effective complex Bose-Hubbard model 
can also emerge even if a system is invariant under gen- 
eralized time reversal transformation T. 

For the sake of simplicity, let us consider a 2D bosonic 
gas in a triangular optical lattice with primitive vectors 
ai and a 2 , see Eq. ((3]), subjected to a periodic driving 



that results in a motion of the lattice on a periodic tra- 
jectory with a velocity 



R(t) = e x cos(ujt) —e v sin(<jt)+ 

muj mio 



2mw v y ' 



(21) 



The Hamiltonian ©, where F(t) = — mR(i), does not 
change under the generalized time-reversal transforma- 
tion T = P X T, where T is the time-reversal operation 
and P x stands for x — > —x. [2^. The Wannier basis vec- 
tors are not invariant under the P X T transformation and 
consequently the Floquet Hamiltonian is not necessar- 
ily represented by a real symmetric matrix in that basis. 
Thus, we may expect that the effective Bose-Hubbard 
Hamiltonian (JTUJ), which is a single block of the entire 



Floquet Hamiltonian matrix |14l . |15| , may possess com- 
plex tunneling amplitudes. Indeed, the effective tunnel- 
ing amplitude 



•h = — dt cxp — ai • R(t) 



T 



J Jo 



F x a 

hijj 



(22) 



is real but 
Ji = 



J 
T 



1 ffl 

dt exp I — — a 2 • R(i) 
o V h 



+ 00 



= J J * 



n— — oo 



(^l_\ p -i2n^ (23) 



hL) Jn \2hw) e 



is complex valued if F y ^ and Q = F x j (\/3Fy) is not 

a multiple of tt/2. In Eq. ^ K x = a^F^ + 3F^/2 

and K2 = a\/3Fy/2. Amplitudes related to tunneling in 
different directions are indicated in Fig. [5] This figure 
shows also that the generalized time-reversal symmetry 
is preserved by the effective Hamiltonian, i.e. complex 
conjugation combined with reflection x — > —x does not 
change the Hamiltonian. 

Now let us consider a problem if the absence of any gen- 
eralized time-reversal invariance of a periodically driven 
lattice system is a sufficient condition in order to obtain 
the effective Bose-Hubbard Hamiltonian with complex 
tunneling amplitudes? We will see that it is not, i.e. 
we are able to construct an example where any general- 
ized time-reversal symmetry is broken but tunneling am- 
plitudes, in an effective Bose-Hubbard Hamiltonian, are 
real valued. Consider a simple ID model, driven with a 
periodic piecewise constant velocity 

v A , < t < T x 
v B , T x < t < 2Ti 

vc, 

-V B , ■■■ ' 

-VA, 

-v c , 5T x <t< 6T1 



R(t) 



(24) 
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where Ti = T/6 and v A > v B > v c > 0, see Fig. [9] It 
is easy to check that the resulting tunneling amplitude is 
real valued, 



_ / dt cxp ( —R( t ) 



(mv A \ , (mv B \ 
cos I — — ) + cos I — - — J 



V h J 



\ h 



(25) 



Effective Bose-Hubbard Hamiltonian is a single block of 
the entire Floquet Hamiltonian of a shaken lattice system 
[Til [H| . We see that such a block can be real even if the 
entire matrix is expected to be complex Hermitian. 

In Ref. [TtJ two conditions have been considered in or- 
der to achieve an effective complex Bose-Hubbard model. 
That is, a tunneling amplitude between two lattice sites 
Yi and Yj can be complex, if the projection of an inertial 
force fij(t) = (Yi — Yj) ■ F(t) breaks: (a) reflection sym- 
metry for suitable time r, i.e. t) = fij{—t—T) and 
(b) shift anti-symmetry, i.e. fij(t — ir/uj) = —fij(t). The 
model presented in Fig. [5] consitutes an example where 
the conditions are broken but the tunneling amplitudes 
remain real valued. 



VI. CONCLUSIONS 

In this paper we propose a complete realization of a 
quantum simulator of the classical XY model in ultra- 
cold atoms in a three-dimensional triangular optical lat- 
tice and present its predictions. In order to manipulate 
independently values of three different couplings we con- 
sider a motion of the lattice constrained to an ellipse. 

Out of simulations we obtain a wide variety of spin 
configurations, some of which being generalizations of 
the two-dimensional triangular lattice model, and a full 
phase diagram. In particular we observe the emergence 



of a new phase, enabling a second order phase transition 
between ferromagnetic and rhomboidal phases, and a dis- 
continuous phase transition between two types of spiral 
phases. This discontinuous phase transition is a conse- 
quence of J3 interaction, absent in the two-dimensional 
model. 

Optical lattice experiments start usually with ferro- 
magnetic state and then parameters of the lattice are 
slowly or suddenly changed in order to switch to a differ- 
ent phase. We have applied Bogoliubov approach which 
allows us to analyze stability of stationary mean-field so- 
lutions under a change of system parameters. For suffi- 
ciently strong particle interactions the initial ferromag- 
netic state collapses towards a ground state — for typical 
experimental parameters it takes place during a few mil- 
lisecond. On the other hand if the interactions are too 
weak the initial state loses its stability but it does not 
evolve towards a system ground state. 

In the present paper we mostly concentrate on a Bose- 
Hubbard model with real tunneling amplitudes. How- 
ever, we consider also a problem of the realization of com- 
plex tunneling amplitudes by means of periodic shaking 
of an optical lattice. It is known that a time periodic per- 
turbation that the breaks time-reversal symmetry leads 
to a Floquet Hamiltonian that is represented by a com- 
plex Hcrmitian matrix and the effective Bose-Hubbard 
model can be expected to possess complex tunneling am- 
plitudes. However, we show that it is also possible to 
realize such amplitudes with a perturbation that is in- 
variant under a time-reversal transformation. 
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